Drifting Spatial Structures in a System with Oppositely Driven Species 
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A system consisting of two conservative, oppositely driven species of particles with excluded vol- 
ume interaction alone is studied on a torus. The system undergoes a phase transition between a 
homogeneous and an inhomogeneous phase, as the particle densities are varied. Focusing on the 
inhomogeneous phase with generally unequal numbers of the two species, the spatial structure is 
found to drift counter-intuitively against the majority species at a constant velocity that depends 
on the external field, system size, and particle densities. Such dependences are derived from a 
coarse-grained continuum theory, and a microscopic mechanism for the drift is explained. With 
virtually no tuning parameter, various theoretical predictions, notably a field-system-size scaling, 
agree extremely well with the simulations. 
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I. INTRODUCTION 

In the recent decade, there was considerable interest 
in the statistical mechanics of a variety of systems in 
stationary, but non-equilibrium, states. Notable exam- 
ples include fast ionic conductors, surface growth, elec- 
tromigration, flux creep in superconductors, propagation 
of defects and cracks, electrophoresis, and granular as 
well as traffic flow. Apart from practical applications, 
the interest lies in the need to establish a sound foun- 
dation for non-equilibrium statistical mechanics, on par 
with the Boltzmann/Gibbs formulation for systems in 
equilibrium. To pursue these goals, many authors have 
proposed simple models, just as Lenz and Ising did y] 
in order to understand the phenomena of phase transi- 
tions of a magnet in thermal equilibrium. Along these 
lines, Katz, Lebowitz and Spohn [Ej introduced the sim- 
ple driven Ising lattice gas, as an "entry" into the physics 
of non-equilibrium steady states. Since then, this field 
has steadily grown, so that there now exists many varia- 
tions and generalizations of the proto model . 

One of the most natural generalizations are systems 
with a second species of particles. The simplest of these is 
a model with equal numbers of oppositely "charged" par- 
ticles, driven by a uniform external "electric" field and 
diffusing on a periodic, square lattice Q. With no in- 
terparticle interactions, except the excluded volume con- 
straint, this system exhibits a phase transition, for criti- 
cal values of the particle density and external field, from 
a homogeneous disordered state with sizable current to 
an inhomogeneous state with minute current. Particles 
of the opposite charge impede each other and "lock up" 
into a dense region. By symmetry, the average location 
of this region is time- independent. Since its inception, 
a number of its properties are reasonably well under- 
stood &-H| , while a variety of related ones are proposed 



|p[-pl|. However, none of these studies focused on a sys- 
tem with unequal numbers of the two species. Assuming 
that a locked- up state still exists, one should not expect 
the dense region to remain stationary. In particular, we 
can expect a larger number of the majority species to 
lie within this region, and so, naturally expect the block 
to advance with the majority. In this paper, we study 
such inhomogeneous states with both simulation and an- 
alytic techniques. Perhaps most surprising of the results 
is that, in the ordered phase, the spatial structure, as a 
whole, drifts in a direction opposite to the intuitive pic- 
ture above! 

The remainder of this paper is organized as follows. 
In the next section, we will provide specifications of this 
model and some details of the Monte Carlo runs. We 
present the simulation results of the counter-intuitive mo- 
tion of inhomogeneous states and suggest its microscopic 
mechanism in Section III. Section IV will be devoted to 
the continuum mean-field approach, which was relatively 
successful in describing the charge-neutral model [0,pp] 
and will be re-analyzed for the more general case here. 
These theoretical predictions are then compared to the 
simulation data in Section V. Particular attention will be 
paid to the scaling of the drive with the system size, and 
the dependence of the drift velocity on control parame- 
ters. We end with some concluding remarks in the last 
Section. 



II. A MODEL FOR DIFFUSION OF TWO, 
OPPOSITELY BIASED, SPECIES 

Generalizing the work of Ising, Potts |lj] and Blume, 
et. al. U^ introduced models which consists of only three 
or more states per site in order to describe various sys- 
tems such as magnets with spin one or higher and ternary 



mixtures. Along similar lines, the natural generalization 
of the driven Ising lattice gas ||] would be models of sev- 
eral species of particles, driven far from equilibrium by 
some "external" field. Clearly, there are many physical 
systems for which such models may be applicable. Here, 
we will focus on the simplest one [Q . 

On a square lattice with Lx x Ly sites, we place N± 
particles with "charge" ±1. At each site, there will be at 
most one particle, regardless of its charge. Thus, a con- 
figuration of our model is completely specified by the set 
of occupation numbers {n± {x,y)}, where n± (x,y) = 1 
or 0, if there is a ± particle at site (a;,y) or not. Apart 
from this excluded volume constraint, there is no interac- 
tion between the particles. However, there is an external 
"electric" field, E, chosen to point in the ~\-y direction, 
so that a -|-(— ) particle is biased against moving in the 
— (+)y direction. Specifically, the system evolves by ran- 
dom updating. In each trial, a pair of nearest neigh- 
bors is randomly chosen. If it is a particle-hole pair, 
then the particle hops into the hole with a probability 
min{l, 6=*=^^"} for the ± species, where a denotes the 
direction of hopping. Lj. x Ly such trials constitute one 
time step (or one sweep). Finally, we impose periodic 
boundary conditions, so that our lattice is in fact a torus. 
For later convenience, let us define the terms "overall 
mass density" and "overall charge density" , given respec- 
tively by 






and 



N. 



N^ 



LirpL/.., 



(1) 



Clearly, for E — Q, there is in fact no distinction be- 
tween the two species. The system is purely diffusive and 
uninteresting. On the other hand, for _£ > 0, particles 
of the opposite charge impede each other. This mutual 
blocking is so severe that the system displays drastically 
different characteristics if the particle densities are high 
enough. In all the previous studies H-I^l of this model, 
q is restricted to zero for simplicity, so that there are 
only two control parameters, {E,m), besides the system 
size. There, for fixed E, say, the steady state of the sys- 
tem is disordered and homogeneous, provided m is small 
enough. By symmetry, the two opposing particle currents 
are the same, on the average. Thus, the (average) hole 
current, C, is zero, while the (average) charge current, J, 
is non-trivial. As m increases, J increases sublinearly, as 
a result of the excluded volume constraint as well as the 
mutual blocking. Once m rises beyond a critical value, 
iriciE), a phase transition occurs so that the system is 
ordered into an inhomogeneous state. In this state, three 
regions can be roughly identified: one particle-poor and 
two particle-rich zones, one of each species. As might be 
expected, these regions span the transverse dimension of 
the lattice (Lx), with each particle-rich zone impeding 
the "forward motion" of the other species. For systems 
with 0(1) aspect ratios, these zones are purely trans- 
verse to the drive, i.e., the densities are homogeneous 
in X. The current drops to vanishing values. If E is 



sufficiently large, this transition is extremely sharp and 
dramatic [Q. With larger aspect ratios, the system often 
locks up into somewhat different states, with zones span- 
ning both X and y, i.e., wrapping around the torus with 
non-trivial winding numbers H . The current still suffers 
a drop, though not to vanishingly small amounts. In ei- 
ther case, once lock-up occurs, these zones are stationary 
on the average, since C = always. 

In this paper, we will study systems with unequal imm- 
bers of the two species. With q ^ 0, many of the previous 
properties will be different, although we still expect the 
presence of a phase transition. For example, in the ho- 
mogeneous state, C will not vanish and propagating fluc- 
tuations are possible. Defering a comprehensive study of 
this model to a later publication |1J], we will focus here 
only on the inhomogeneous state, in which the zones are 
expected to drift. As we will demonstrate and explain, 
the system displays a counter- intuitive feature, i.e., the 
spatial structures drift in the direction favored by the 
minority species. For example, the inhomogeneities will 
drift in the negative y direction if q is positive! 



III. DRIFTING STEADY-STATE STRUCTURES 
FROM SIMULATIONS 

Since our purpose here is to study the novel properties 
associated with unequal numbers of the two species, we 
carry out simulations with fixed particle density (p+ = 
N^/ LxLy) for the + species and varying density for the 
— species (p- = N-jLxLy < p+), corresponding to 

q = P+- P- > ■ 

The behavior for the case of g < may be deduced 
simply by symmetry. Specifically, we choose p+ — 1/4, 
Lx = 10,20,40, Ly = 40,160,320, and E ranging from 
about 0.1 to 1. These parameters are chosen in order 
to probe the q > region close to the q — inhomo- 
geneous states near the transition mentioned above, for 
we expect the novel properties to be more pronounced 
there. In contrast, the particles can hardly move deep 
in the locked-up phase at higher densities. Simulations 
for different Lx^s show that the effect of Lx is negligible, 
as found for the symmetric, q = case 0. Thus, unless 
La: ^ Lj, |q| which we will not consider here, we are deal- 
ing with a system in which only one of the dimensions 
plays an essential role. 

Starting from the inhomogeneous state with p_ = p+, 
we find a phase transition into a homogeneous state as we 
gradually decrease p_ with p_|_ held fixed. On the (g, m) 
plane, the phase boundary between the homogeneous and 
inhomogeneous states may be located this way, which is 
symmetric about the m axis. However, a detailed dis- 
cussion of the phase diagram is beyond the scope of this 
paper. 

Focusing on the properties of the q > inhomoge- 
neous states, we find that the locked-up region of the 



two species drifts backwards with respect to the driving 
direction for the majority (+) species at a definite ve- 
focity V that increases with q but decreases with E, as 
shown in Fig. |l|. Fig. g shows a typical inhomogeneous 
configuration in steady state for g > 0. The steady-state 
ensemble averages of the local density profiles for the two 
species, p-(-(y,t) and p_(?/,i), are measured. Due to the 
drift, they are functions oi u ^ y — vt alone. Fig. ^ dis- 
plays the steady-state density profiles p±{u) for various 
values of g > 0. 

To understand the microscopic mechanism for this 
backward motion, it is instructive to consider the role 
of holes inside the two particle-rich zones. The probabil- 
ity for a hole to diffuse against the drive into these zones 
from the outer zone boundaries is suppressed by E via 
e~^. Thus, provided E is not too large, there are finite 
densities of holes inside the block. For q — Q, these den- 
sities are on the average the same in the -|- and — zone. 
For g > 0, there are more holes in the — zone because it 
is thinner (its thickness given roughly by p-Ly). With 
holes available on the inner (-1 — ) zone boundary, a parti- 
cle may escape from the block to the particle-poor region 
through the zone of the opposite species. Driven along E 
then, it eventually returns to the outer edge of its own 
zone due to periodic boundary conditions. When a par- 
ticle leaves the inner zone boundary, a hole is left behind 
which may drift in cither direction towards an outer zone 
boundary, returning to the hole-rich region. For q = 0, on 
the average, the number of holes impinging on an outer 
zone boundary equals the number of incoming particles. 
Thus, apart from a migration of particles from the inner 
to the outer zone edges, the cluster remains stationary. 
For g > 0, however, it is relatively easier for + particles 
to migrate. It results in more particles than holes im- 
pinging on the outer -I- zone boundary, and the opposite 
for the — zone. It is this imbalance that causes the whole 
cluster to drift backwards. Of course it is clear from our 
argument that v must vanish ii E — oo. 

In order to see how this arises theoretically and to ex- 
plore the E and q dependence of u, we now turn our 
attention to a continuum description. 



IV. A CONTINUUM MEAN-FIELD 
DESCRIPTION 

Following previous studies of this model |0~||] , we rely 
on a mean-field type continuum theory to understand 
the macroscopic properties here. The equations of mo- 
tion for the densities, first proposed in Ref. Q, need no 
modification and, for completeness, summarized below 
(Section III. A). However, the overall constraint on the 
charge density will be different, leading to qualitatively 
new behavior such as drifting inhomogeneous solutions 
(Section III.B). 



A. Equations of Motion 

To describe the long- wavelength, low-frequency behav- 
ior of our model, we make use of the continuum approach, 
in which the discrete variables of both the lattice and the 
occupation numbers, n± (a;, y) , are replaced by continu- 
ous ones for the densities and space-time: p± (r,i). For 
r, we will continue to write {x,y), which should not lead 
to any confusion, and let x e [0, L^), etc. The evolu- 
tion equations of these densities may be "derived" by 
taking the continuum limits of the mean-field approxi- 
mation to the Master equation |IQ|; or they may simply 
be postulated through considerations of symmetries. In 
the former approach, the parameters in the continuum 
equation can be related to the microscopic rates. Since 
we will not be concerned with the absolute time scale, 
one parameter may be absorbed into the definition of t. 
In other words, we will set the diffusion constant, for the 
undriven case, to be unity. Only one parameter remains, 
associated with the driving field. If the naive continuum 
limit approach is taken, then it is 



8 = 2 tanh(S/2) 



(2) 



With these considerations, we study the following 
equations of motion, written in the form of continuity 
equations: 



dp± 
dt 



V • [</>Vp± - p±V(/) T p±0f y] 



(3) 



where (j) = 1 — /?+ — p_ is the density of holes. Notice 
that the first two terms in these equations describe free 
diffusion of two distinguishable species of particles. The 
last term corresponds to the Ohmic currents, with /9±0 
being the usual density dependent conductivity. It is 
also natural to consider the sum and difference of these 
equations. Defining ip = p^ — p_ to be the charge density, 
they take the form 






(4) 

(5) 



These are precisely the equations in Ref. [0 . To apply 
to our problem, we only need to impose 



LxLy 



ip dxdy = g > 



instead of g = 0. The other constraint, 
1 



LxLy 



dxdy = 1 — TO , 



(6) 



(7) 



as well as the periodic boundary conditions for the den- 
sities, of course remain unchanged. 

We may simplify these equations further, by absorbing 
£ into the scale of y. There is no need to write new 



equations, since we can simply drop £ from (y,^ while 
keeping in mind that Ly must be replaced by £Ly. That 
the drive provides an intrinsic length scale implies that 
the ordinary thermodynamic limit {Ly — > oo) must be 
taken along with £ — + while holding the product £Ly 
fixed. However, this simplification may be too confusing 
and will not be used here. Due to the central role played 
by £Ly, let us define 

e ^ £Ly (8) 

for future convenience. These equations, (Q-0), com- 
pletely specify the dynamics of our model. 



B. Inhomogeneous Steady States 

Next, we study steady state solutions to these equa- 
tions which are spatially mhomogeneous. Since we ex- 
pect the densities to be time dependent, let us seek so- 
lutions with a constant velocity, w, namely, 4'{x,y — vt) 
and ip{x,y — vt). Simplifying further, we note that all 
the states we observed in simulations are homogeneous 
in X, so that we will restrict ourselves to functions of the 
form: 

(j){u) and -0(u) , (9) 

where u = y — vt. Inserting these into (HH), we have 

—vd(j> = d [d(j> + 4>i}£] and 

—vdijj = d [(j)dip — tpdcj) — (/)(1 — <j))£] , 

where d stands for d/du. Integrating once, we obtain 



d(j) + £<j)'i{j+v(j) = ~C 
(pd^ - ^pd(j} - £(j){l - (j))+vi) = -J 



(10) 



The constants, C and J, may be interpreted as the two 
steady state currents for the holes and the charges respec- 
tively in the moving frame. In the "lab-frame" , the cur- 
rents should be mhomogeneous, due to the anticipated 
drift of the block. As it stands, equations (llG) contain 
three unknown constants (w, C, J), which will have to be 
fixed by three conditions, i.e., solutions be of period Lj,, 
(^ and (0). However, analytically, {v,C,J) appear to 
play more the role of control parameters while (Lj,, g, m) 
will emerge at the end. In this way, the analytic approach 
is somewhat opposite to that of simulations, where the 
latter (former) are the control (dependent) variables. 

To find the solutions, we follow previous studies and 
introduce variables which simplify the structure of these 
equations: 



X = l/i/i and ijj = ipx 



(11) 



Note that, unlike the physical densities which are 
bounded, x ^ [IjOo] and IV'I < X — 1 G [0, oo] . Now, 
( p^ ) becomes 



dx = £'4' +VX + Cx^ 

d^ = ^(X - 1) - Jx^ - viiJx 



(12) 
(13) 



Eliminating ^p, we again arrive at an ordinary differential 
equation for only one variable: 



X 



(x - 1) + Jx = V X + vCx 



"X 1 - 



Cx' 



(14) 



where J = J/£, etc. Also, prime denotes d/ du£, showing 
again the central role played by £ in setting the length 
scale. For clarity, we have placed on the right hand side 
of this equation all the extra terms due to q ^ 0. 

Unlike in the neutral system, the interpretation of (|lj) 
as a particle "moving" in a potential has to be modi- 
fied, since there are "velocity" (i.e., x') dependent "force" 
terms. In general, periodic "motion" would be impossi- 
ble. Of course, here, we must insist on the existence of 
such solutions. The consequence is a constraint on the 
last term in (|l4|). In particular, multiplying this equation 
by x' and integrating over the full period, we are led to 

/(x')^ u (1 — x) + SCx du =■ 0. Since we are concerned 

with inhomogeneous states, we can expect (x')^ to be 
positive, except for isolated points. Thus, it is possible 
to interpret {x')'^d.u/ J{x')'^du as a new measure on the 
interval u G [0, Ly] and define a new type of average: 



(•) 



jix'rdu 



(15) 



Using this notation, we may write a simple relationship 
between the drift velocity and the hole current in the 
moving frame 



v{x-l)=2C (x) 



(16) 



Since x > 1 typically, we conclude that the drift is in the 
same direction as the hole current. A similar relationship 
can be found by integrating Eq. (|lj), after multiplication 
by (j). The result is: 

-q£ =v + Cx , (17) 

where bar is the normal average: 

f • du 
• = . 

Ly 

Eliminating C between (nm and (0), we see that 

, (x-i)x 



(18) 



V 



2(X) 



(19) 



is negative definite. So, for example, if the majority 
species is positive (i.e., more particles are driven "up- 
wards"), then the drift of a block state will be "down- 
wards" ! This behavior is, naively, quite surprising, since 



we expect the particle-rich zone to contain more particles 
of the majority species so that the entire block should 
"advance" with the majority. Instead, the block drifts 
in the opposite direction. On closer examination, we find 
that, since the negative region (in this example for q > 0) 
is thinner, it is easier for positive particles to get through 
the blockage. Then, due to the periodic boundary condi- 
tions, these particles pile up "behind" the positive region. 
As a result, the entire block appears to drift "backwards." 
This picture simply provides another perspective on the 
intutive arguements in Section III. Here, we have proved 
that the structure "retrogrades." 

Clearly, this analysis also leads to C/q being negative, 
i.e., the holes moving contrary to the majority species, 
which is hardly surprising. Before closing this subsection, 
we should comment on a number of other constraints on 
the three unknowns (w, C, J), independent of the specific 
values of (Ly, q, m). 

In order to have any periodic solution at all, there 
must be some form of restoring "force" in (O) . Exam- 
ining the "potential" part of this "force", i.e., (x ~ 1) ~ 

J — v^j x^+^C'x'^, we see that there would be no "well" 



to trap the particle, unless 

J>v^ . (20) 

On the other hand, at least J — v^ < 1/4 is needed, even 
in the neutral case. The cubic term further exacerbates 
the situation. The constraint that a "well" exists turns 
out to be 



DC 



4-18 J- 



[J-,') 



27vC 



< J - v^ 



(21) 



More information can be gleaned from regarding JT^ ) 
and (03) as flows in a "phase" plane. If a periodic solu- 
tion exists, it would correspond to a closed loop and, by 
continuity, there would be, generically, at least one focus 
(fixed point with spiral orbits) lying within. In order to 
have a physical solution, this fixed point must lie in the 
physical region: x G [I7O0] and jV'l < X ^ 1- Since any 
fixed point must lie on the curve V' = ~ ( "^X + C'x^ 
obtain 



{l + vf > -4(7 



(22) 



Recalling that both v and C are negative, this confines 
both to be small quantities. 

In our case, it is easy to check that there are three fixed 
points, one of which always lie outside the physical re- 
gion. Of the remaining two, one is a focus and the other, 
a saddle. Based on the characteristics of the neutral sys- 
tem, we expect our solution curve to run in between these 
two fixed points. Unlike the neutral case, however, the 
flow is not Hamiltonian in general and, in particular, the 
focus is not necessarily a center (i.e., eigenvalues corre- 
sponding to the flow linearized about the fixed point not 



necessarily pure imaginary). Thus, there is no guaran- 
tee that we can find a periodic solution. One possible 
scenario is that a unique limit cycle exists for any given 
{v,C, J), provided they respect the inequahties (|20|-P2|). 
Another is that {v, C, J) satisfy a specific relation which 
allows for the existence of periodic solutions. A natural 
constraint to impose is that this fixed point be a cen- 
ter, with its associated eigenvalues being pure imaginary. 
This condition is equivalent to setting the coefficient of 
the last x' term in ( [l4|) to zero at that fixed point. This 
gives us an additional formula for the velocity: 



2C_ 

i_ 

X* 



(23) 



where y* denotes the fixed-point value of x at the center. 
Using (|l^) , it is easy to see that the value of the densities 
at any fixed point satisfies a cubic equation with param- 
eters {v, C, J). Eliminating x* by (|23|), we then obtain a 
quintic algebraic equation for v alone: 



8C' 



32C,, 
3Jv^ 



^v- 



'i2C„y + 2 JC„ 
- 2v^ = 0, 



18u-^ 



UrnV 



(24) 



where Cm = ^v — C is the mass current in the moving 
frame. However, though this scenario guarantees peri- 
odic solutions at the lowest order in the neighborhood of 
the fixed point, we are unable to prove that, beyond the 
linear level, this condition is either necessary or sufficient 
for the existence of periodic solutions. There is, never- 
theless, some numerical evidence that such solutions are 
available, as we will show in the next Section. Equation 
( p4| ) prescribes a surface v{C, J) in the C — J plane. Sub- 
ject to numerical uncertainties, we find that the param- 
eters (f, C, J) generated by simulations indeed span a 
surface consistent with this scenario. 



V. FIELD-SIZE SCALING AND COMPARISONS 
WITH SIMULATIONS 

To determine how close the continuum model corre- 
sponds to the discrete model, we subject our theoretical 
predictions to the tests of simulations. The first is con- 
cerned with the scaling behavior in the system size and 
field strength. Choosing the alternative set of control 
parameters (Ly, q, m) in favor of (u, C, J), Eqs. (|l2|) and 
(O) implies that the solutions for the densities obey a 
simple scaling form (cf. [p||qjl0||): 

x(u, E, Ly, q, ra) = F^{£u, £Ly, q, m) 

= F^{u/Ly,£Ly,q,m), (25) 

where F^ and F^ are appropriate scaling functions. Sim- 
ilar scaling form holds for V', of course. Monte Carlo 
simulations, with fixed p+ and varying p_, using a wide 
range of Ly and E with fixed e = £Ly show excellent 



agreement. An example is shown in Fig. Q. One immedi- 
ate implication of this result is that the thermodynamic 
limit has to be taken with care, as phase transitions sur- 
vive only in the double limits i? — > and Lj^ — > oo with 
£ held fixed. 

A more stringent test is to check to what extent the 
data actually satisfy the differential equations (|lj) and 
(|l3|). In the continuum description, it is more natural to 
use the currents in the moving frame because they are the 
integration constants. In simulations, the (spatial and 
temporal) average currents in the lab frame are more ac- 
cessible. They are related. For example, the hole current 
in the lab frame is given by C -I- vcj), which is greater than 
C, in magnitude. With no tuning parameter, the equa- 
tion (nj) for the hole density fits the data very well (see 
Fig. |5|(a)), but there are appreciable discrepancies in the 
equation ( |l3| ) for the charge density in the particle-rich 
region (see dashed line in Fig. H(b)). Similar discrepan- 
cies were also observed in a closely related model con- 
sisting of two species driven along orthogonal directions 
[Q. In that model, the asymmetry between the (H — ) 
and ( — h) nearest-neighbor correlations along the field di- 
rection was shown to give rise to additional cubic terms 
of the form ±Xp^p-(f>£y in the currents for the ± species, 
which enter inside the brackets on the RHS of Eq. (||). 
Due to the opposite signs, they cancel out in the equation 
(0) for the hole (or mass) density but contribute an extra 
term to equation (||) for the charge density. These terms 
represent the lowest order corrections to our mean-field 
equations in Section IV. After such a term —2X£p^p-x is 
added to the RHS of Eq. (O) , significant improvements 
are found, for a suitable proportional constant A (solid 
line in Fig. B(b)). 

Further comparisons are concerned with the mean cur- 
rent and drift velocity. The mass current —C — v(f> in the 
lab frame, finite for q > 0, is simply given hy MS, which 
is obtained by integrating the first of Eq. (|lO|). Excel- 
lent agreement with simulations is found, as shown in 
Fig. |6|(a). This comparison does not involve the A cor- 
rection terms. Other predictions, such as (|l9|) and (23), 
derived by using both the hole and charge equations with- 
out A, do not agree as well. Fig. ^(b) exhibits increasing 
deviations as q increases. The agreement, however, can 
again be significantly improved by including the A terms, 
with the same choice of A w 1.5 as in Fig. ^. With A 7^ 0, 
(p^) is slightly modified: 



^-'^ 



2X+{x-l)^X 
2(1 - A) (x);, 



where, similar to (p^). 



lix'/x^rdu 



Also, 



becomes 



(2-A)C 



1 



A-l 

X* 



(26) 



(27) 



(28) 



where y* is approximated by the spatial maximum Xmax 
in Fig. |g, the error incurred is very small as both quan- 
tities are much greater than 1. Detailed discussions and 
derivations of the A terms will be presented elsewhere 

The final convincing evidence for the quantitative 
agreement is a direct comparison of the density profiles. 
For simplicity, we consider only the case of A = 0. The 
fixed-point condition mentioned near the end of Sec. IV 
picks out a unique v for a given set of currents (J, C) 
via (p5). Equations ( p^ ) then contain no free parameter, 
and we can obtain the profiles by numerical integration, 
using for instance the Runge-Kutta method jlsj . A typi- 
cal comparison with the simulated profiles using the same 
set of parameters (g, m, Ly, E) is presented in Fig. |7|. The 
agreement for this case of rather small q w 0.09 is again 
excellent, even without A. We expect more deviations for 
larger q, where the A terms can no longer be ignored. 

These comparisons provide strong support to our claim 
that the continuum model, which may be systematically 
refined if necessary, represents a surprisingly accurate de- 
scription of the simulated discrete model. 



VI. CONCLUDING REMARKS 

To summarize, we have studied, using both Monte 
Carlo techniques and the continuum meanfield method, 
a diffusive system of two species of particles, driven in 
opposite directions by an external field. With periodic 
boundary conditions imposed, this system settles into a 
non- equilibrium state with a steady current. For simplic- 
ity, we have restricted ourselves to non-interacting parti- 
cles, apart from an excluded volume contraint. Thus, as 
the particle densities increase, this system undergoes a 
phase transition, from a homogeneous disordered phase 
with a high current to an inhomogeneous one with minute 
current. In the latter state, the two species impede each 
other so much that they form a blockage of high, local 
particle density. When the particle densities are unequal, 
this spatial structure displays a counter-intuitive behav- 
ior. It drifts with a constant velocity, in a direction op- 
posite to that favored by the majority species. Remark- 
ably, simulation results agree reasonably well with most 
aspects of the theory, especially the prediction of £Ly 
scaling. 

On the other hand, there are clear signs of disagree- 
ment, mainly in connection with dl3). Since our the- 
ory is based on meanfield assumptions, one avenue for 
improvement is to take some correlations into account. 
The simplest addition, involving cubic terms prt in (H) 
lead to significant improvements. Encouraged by these 
findings, we are undertaking a comprehensive study, in- 
cluding a general phase diagram in the {q, m) plane, of 
the effects of such terms. In this paper, we have focused 
only on the drifting inhomogeneous state. Although we 
have performed some analysis for the homogeneous state 



p4[ , much remains to be investigated. For example, it 
would be desirable to observe, in simulations, the drift of 
fluctuations from the uniform densities. Of course, as in 
the neutral case, we should expect long range correlations 

When restricted to one-dimension (i.e., one column), 
this model is exactly soluable ||lj,|l^, since the order of 
any particular string of -l-'s and — 's is invariant and no 
phase transitions can occur. With open boundary condi- 
tions, it can be mapped onto the Rubinstein-Duke model 
for electrophoresis p8| , and more interesting phenomena 
can occur. Beyond the simple model studied here, there 
are many other generalizations which may be relevant to 
a variety of physical systems. We mention only a few 
here. 

The existence of the inhomogeneous state depends cru- 
cially on the mutual blocking between the species. To 
find out the importance of this effect, we may introduce 
"charge exchange" processes which takes place at a frac- 
tion of the particle-hole exchange rate. As in the neutral 
case [Q, we can expect to find the transition between 
the disordered and the inhomogeneous states to be both 
continuous and discontinuous. It would be interesting 
to map out a complete phase diagram in the {£Ly, m, q) 
space. Further, such a system can display interesting be- 
havior even in one dimension ||19[, especially with open 
boundaries pG]. In particular, it is closely related to the 
model of "first and second class particles" (e.g., cars and 
trucks) moving on a ring, in which a first class parti- 
cle is allowed to overtake a second class one with some 
rate uM. In that case, there are distinct phases, with 
properties reminiscent of our inhomogeneous states. The 
qualitative vs. quantitative similarities should be inves- 
tigated. 

Another generalization involves the two species being 
driven in orthogonal directions. These models are moti- 
vated by the phenomena of traffic flow in c ity blocks and 
display a considerable variety of phases [p|jlC|]. However, 
we believe that there are no studies with unequal num- 
bers of the two species (though we are aware of a study 
with varying densities of a third species [^). As we have 
shown in this paper, it is likely that novel behavior will 
be found if the species are not exactly balanced, which 
in view of the physical motivations should be the more 
generic cases to study. 

In all the models mentioned so far, there is no inter- 
particle interaction, except for the excluded volume con- 
straint. It is clearly important to ask what the effects 
of including such interactions are. In particular, even in 
the absence of external drives, there would be rich phase 
diagrams if interactions were present [ [isHl^ . Thus, it 
is natural to inquire how the drive would modify these 
phase transitions. We are aware of only one study of a 
driven system with two interacting species |23[| . Though 
the regime investigated was extremely limited, several 
novel features were found. 

Finally, we point out that, in physical systems such as 
fast ionic conductors, the two species may be of different 



mobilities and different "charges" . These properties add 
two entirely new dimensions to the phase space, leading 
to seemingly endless horizons for future explorations. 
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FIG. 1. The position of the center of mass j/cm vs time for 
the inhomogeneous state. The steady, backward drift veloc- 
ity decreases with increasing average density of the minority 
phase. p+ — 0.25 and E = 0.1976 except otherwise stated. 
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FIG. 3. Steady-state density profiles of a 20 x 160 system 
for difi'erent p_, at fixed p+ — 1/4 and E — 0.1976. p+{u) is 
on the left, p-(u) on the right, and u — y — vt. 

Leung-Zia PRE-Fig4 




+ E=0.8327, Ly=40 
o E=0.1976, Ly=160. 
- E=0.0986, L =320 




0.0 0.2 0.4 0.6 
U/L. 



FIG. 4. Predicted field-size scaling is confirmed by sim- 
ulations for different E and Ly with fixed £Ly. Lx ~ 20, 
p+ = 1/4, p- = 0.16. 



FIG. 2. A typical 40 x 160 configuration in steady 

state showing the blockage between the two species. The 
open/filled squares represent the upward/downward, or +/— 
drifting species. Note that there are more + particles es- 
caping through the blockage to cause the structure to drift 
downwards. Here p+ — 1/4, p_ — 0.1, E = 0.1976. 
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FIG. 5. Typical tests of local properties of the continuum 
model against simulations for (a) the hole equation (12), and 
(b) the charge equation (13), demonstrating the significance 
of the correction term oc \£p+p-X in the latter. L^ ~ 40, 
Ly = 160, E = 0.1976, p+ = 1/4, and p_ = 0.1. 
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FIG. 7. Excellent agreement between theory and sim- 



ulation for the density profiles. Ly 
p+ = 0.272, p- = 0.184 and A = 0. 
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FIG. 6. Comparison between theory and simulation for 
(a) the mass current in the lab frame, —C — wj), and (b) the 
(negative) drift velocity with and without the correction term, 
as a function of the average density of the minority species. 
La; = 20, Ly = 160, E = 0.1976, and p+ = 1/4. 



